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On  the  Formulation  and  Numerical  Evaluation  of  a  Set  of 
Two-Phase  Flow  Equations  Modelling  the  Cooldown  Process 

Stephen  Jarvis,  Jr. 

A  model  of  transient  two -phase  flow  in  a  pipe  is  con- 
structed in  Eulerian  coordinates  assuming  a  single  velocity 
but  independent  temperatures  in  the  two  phases .  Experiments 
on  the  numerical  integration  of  the  system  for  Cooldown 
problem  by  both  lax  and  Courant-Isaacson-Rees  methods  indi- 
cate that  a  very  fine  spatial  difference  net  must  be  used  to 
compensate  for  the  numerical  diffusion  essential  to  computa- 
tional stability  if  a  second  surge  is  to  be  realized. 

1.   Introduction 

This  report  is  intended  to  summarize  progress  to  date  on  the  formu- 
lation and  numerical  integration  of  systems  of  equations  modelling  general 
two-phase  flows  in  a  pipe.  Principal  application  would  be  to  the  study 
of  the  development  of  choked  flows  in  tubes  of  varying  cross -section,  and 
to  the  Cooldown  problem.  The  latter  problem,  the  cooling  of  warm  transfer 
lines  resulting  from  the  flow  of  cryogenic  liquids,  has  been  the  target 
of  the  present  study. 

In  a  private  communication,  a  set  of  Cooldown  equations  was  proposed 
based  on  conservation  laws  for  a  single  fluid,  and  three  distinct  phases 
were  admitted:  all  gas,  all  (compressible)  liquid,  and  a  two -phase  medium 
based  on  a  quality  weighted  average  of  gas  and  liquid  at  a  single  velocity 
(V),  pressure  (p),  and  temperature  (t) .  State  changes  in  the  two-phase 
domain  were  confined  to  the  vapor -pressure  curve  with  quality  (y)  satis- 
fying 0  <  y  <  1.  The  complicated  thermodynamic  structure  of  the  fluid 
(N_)  was  approximated  by  least  squares  fit  to  empirical  curves.  The 
numerical  approach  taken,  Picone's  method,  involved  polynomial  approx- 
imation in  x  and  t  using  the  governing  equations  to  determine  coefficients. 
This  method  is  not  well  suited,  however,  to  hyperbolic  systems  of  partial 
differential  equations  (PDE^)  whose  solutions  may  have  discontinuous 
derivatives  at  interior  points  (and  in  the  case  of  non-linear  equations, 
even  discontinuous  variables,  or  shocks),  although  it  is  applicable  to 
elliptic  and  parabolic  systems  for  which  a  high  degree  of  smoothness  in 
the  solutions  is  the  rule .  Attempts  to  force  polynomial  fits  to 


discontinuities  results  in  highly  unstable  "behavior,  which  was,  indeed, 
observed.  Finite  difference  schemes  are  generally  more  successful  with 
hyperbolic  systems . 

In  addition  to  the  restrictive  assumption  of  homogeneity  on  the  two- 
phase  medium  in  the  above  model  the  location  of  phase  boundaries  is  a 
serious  problem.   If  Ax  is  the  spatial  difference  interval,  it  is  useless 
to  attempt  to  resolve  phase  domains  on  a  scale  even  as  small  as  Ax.  A 
reasonable  value  of  Ax,  based  on  computing  time,  is  10  feet  for  the 
Cooldown  problem,  so  that  it  is  clear  that  if  the  flow  is  not  insensitive 
to  the  structure'  of  phase  domains  of  this  scale,  the  problem  is  compu- 
tationally impossible.   Clearly  the  small  scale  spontaneous  generation 
and  degeneration  of  two -phase  pockets  is  beyond  the  scope  of  the  numerical 
approach  in  general,  and  broad  scale  averages  are  at  all  times  under 
consideration . 

A  fundamental  property  of  continuous  hyperbolic  systems  is  that  the 
flow  at  a  point  (x,t)  is  completely  determined  by  the  flow  at  time  t  -  At 
on  the  interval  (x  -  XAt,  x  +  XAt)  (where  At  is  small),  the  so-called 
"domain  of  dependence".  Here  X  is  the  highest  signal  speed  in  the  fluid 
at  (x,t) .  For  given  Ax  then,  the  smaller  is  X,  the  shorter  is  the 
computation  time  required  to  integrate  over  a  given  real  time  interval. 
Since  liquids  in  general  have  very  high  signal  speeds,  it  is  obvious 
that  compressibility  will  be  a  major  factor  in  the  computing  economics . 
Furthermore,  few  two-phase  problems  of  interests  will  involve  "hammer" 
(liquid-shock),  so  that  considering  the  liquid  to  be  incompressible 
would  be  both  reasonable  and  economical.  A  further  simplification  that 
can  be  made  which  will  probably  not  be  more  significant  than  those  already 
discussed  is  to  assume  the  gas  to  be  a  perfect  gas  and  approximate  the 
vapor -pressure  curve  by  a  straight  line  to  reduce  computation. 

The  flow  of  incompressible  liquid  (having  infinite  signal  velocity) 
can  only  be  handled  in  the  finite  difference  scheme  by  considering  the 
dynamics  of  complete  liquid  columns  bounded  at  each  end  by  compressible 
media.   Thus  the  problem  of  scale  length  and  phase  boundaries  again  appears 


Since  we  are  in  any  case  restricted  to  the  consideration  of  two- 
phase  flows  which  do  not  depend  critically  on  "phase  resolution, "  it 
appears  that  the  whole  matter  could  "be  greatly  simplified  by  allowing 
only  two -phase  flow  everywhere,  0<  Ym^Y  ^  Y*   <  /  ,  placing  then  a 
minimum  and  a  maximum  on  the  quality  y(x,t) .  A  small  volume  of  gas 
then,  while  storing  insignificant  quantities  of  energy  and  momentum, 
adds  compressibility  to  the  "all-liquid"  flow,  Y   -     Yi      •     A.  suffi- 
ciently small  quantity  of  liquid  can  similarly  have  no  effect  on  the 
"all-gas "  flow,  y   -  /u  . 

Taking  this  approach,  we  must  generalize  our  models  now  to  admit 
both  a  liquid  and  a  gas  temperature  to  avoid  being  confined  to  the  vapor 
pressure  curve.  This  extension  appears  to  be  a  worthwhile  generalization 
in  that  heat  transfer  between  gas  and  liquid  is  controllable,  while  super- 
cooling and  super-heating  are  admissible.  A  system  of  five  first-order 
partial  differential  equations  (system  "1V2T")  is  needed  to  describe  the 
flow,  involving  a  single  velocity  V"(x,t),  pressure  p(x,t),  gas  density 
r(x,t),  gas  temperatures  ><,   »  *,  */  ,  and  liquid  temperatures   /^  (  *,  v. 
(Longitudinal  heat  conduction  in  the  wall  is  ignored  also,  so  that  the 
wall  temperature  is  controlled  by  a  single  first -order  equation  in  time 
coupled  to  the  five  mentioned  above  only  through  source  terms  and  is  of 
no  mathematical  consequence) . 

For  further  generality  in  a  one -dimensional  model  a  6th  equation 
was  added  (forming  system  "2V2T")  to  allow  both  a  gas  velocity  V"(x,t) 
and  a  liquid  velocity  U(x,t),  and  requiring  a  description  of  momentum 
coupling  between  the  two  fluids.  But  unless  the  difference  I ^ '  W  is 
adequately  restricted  (by  friction,  for  example),  the  initial  value 
problem  for  this  system  is  not  "well-posed",  i.e.,  there  are  fundamental 
hydrodynamic  stability  problems  involved.   (This  is  evidenced  by  complex 
characteristics).  This  6th  order  system  was  dropped,  then,  in  favor  of 
the  5th  order  (one  velocity)  system  pending  a  study  of  the  meaning  of 
this  instability. 

In  any  fluid  system,  two  natural  coordinate  systems  are  available: 
the  system  fixed  in  space  (fixed  in  the  pipe)  is  called  the  Eulerian 


system,  while  a  system  fixed  in  the  fluid  is  called  Lagrangian.   Only  in 
one -dimensional  transient  flows  does  the  Lagrangian  approach  offer  some 
advantages  over  the  Eulerian,  although  the  fixed  boundary  conditions  (at 
the  pipe  ends)  of  the  Eulerian  scheme  go  over  into  variable  boundary  con- 
ditions .   From  the  coding  point  of  view,  the  finite  difference  method  in 
Lagrangian  coordinates  has  the  great  disadvantage  that  the  difference  net 
shifts  in  and  out  of  the  pipe  and  may  get  too  dense  or  too  sparse  to  admit 
good  approximation  to  the  dynamical  equations.  There  is  involved  frequent 
interpolation  and  renumbering  of  net  points .   Furthermore,  in  the  (2V2T) 
system,  the  occurrence  of  two  velocities  complicates  the  Lagrangian  system, 
(and  might  eliminate  the  natural  advantage  of  the  Lagrangian  system) . 

On  the  other  hand,  the  Lagrangian  schemes  would  have  minimized  the 
difficulty  caused  by  the  great  variations  in  density  and  momentum  density 
that  the  two  phase  system  presents,  and  might  have  described  the  movement 
of  "liquid  points"  without  the  large  mass  diffusion  which  turns  out  to  be 
the  most  serious  problem  of  the  Eulerian  approach. 

Experimentation  with  a  Lagrangian  formulation  would  be  very  desir- 
able.  It  would  involve  no  new  principles,  but  some  tedious  coding 
difficulties . 

2.  Basic  Equations  (2V2T) 

We  shall  assume  the  pipe  has  a  variable  cross-sectional  area  A(x) 
shared  between  a  gas  stream  of  cross -sectional  area  zA  and  an  incompress- 
ible liquid  stream  of  cross -sectional  area  (l-z)A,  so  that  z(x,t)  is  a 
volumetric  quality. 

Let  s  be  the  constant  liquid  density,  and  E  and  F  be  specific  ener- 
gies of  gas  and  liquid. 

In  terms  of  source  strength/unit  length,  S.,  the  conservation  laws 

J 
are  readily  written  for  the  two  streams . 

Liquid  Continuity, 

(A(i-b)s)     +    (A(i-i)s  u)     ~  S  a) 


Gas  Continuity, 


(Air)     f  (AzrV)  x-  Sa 

Liquid  Momentum, 

(A0-*)sU)t  +(A0-i)sU>  +  AO-'Jyfl, 

Gas  Momentum, 

(AirV)jt    *■  (AzrV*  +   AzflK 

Liquid  Energy, 

[A('-*)*(iU\F)] 
+[A0-*UU(iut*F)+A(i-JOf,]     ~  Sr 


Gas  Energy, 


[A*r(J  V3tE)} 

+  [AitV(i  V\E)  +  AzV+]-  S< 


Wall  Temperature, 


(2) 


(3) 


00 


(5) 


(6) 


T    -  S  (7) 

y,  t  l 

We  assume  a  state  equation  for  the  gas, 


Calling  _e^  and  ^c-  the  mass  rates  of  evaporation  and  condensation  per 
unit  length  of  pipe  at  (x,t),  we  have 

o   =  ^  -  *> 

and  (9) 

^  a. 

Let  Jc       £e        Jc  he  perimeters  of  contact  "between  gas -wall,  liquid- 

wall,  and  gas -liquid,  respectively,  at  the  section  (x,t),  and  -?       t^     "X 
the  associated  shear  stresses  (wall  on  gas,  wall  on  liquid,  liquid  on 
gas).  Assume  gravitational  acceleration  Q(x)  in  positive  -  x  direc- 
tion given.   The  scl.     and  -<z,     are  Von  Neuman  -  Richtmyer  non-linear 
pseudo -viscous  forces  introduced  to  permit  shocks,    and  are  described 
later.   Then  we  set: 

S3,  -  A  s  (i-  z)  y  -  r£  2£   +   r0  i. 


(io) 


and 


S,r  &<-*%-  ^£i  -r<£< 


(id 


f  .*  V -  /*  U  -  A  2  -f^t 

Let  Q        Q o     be  heat  transfer  rate  from  wall-into-gas  and  liquid-into- 
gas  per  unit  area,  and  Q     the  wall-liquid  rate.   Then: 


-  j£ 


(12) 


and 


V  *"  VfA  +fA  -  1,-*, 


We  have  taken  the  following  point  of  view  in  formulating  these  sources . 
When  a  quantity  of  liquid  JL  A*-     ±s   evaporated,  its  specific  momentum 
and  specific  enthalpy  are  "brought  to  those  values  in  the  gas  stream 
at  the  expense  of  the  liquid  stream,  and  contrariwise  for  condensation. 
For  example,  the  gas  momentum  source  term  is  -£  V -  ^c  LJ       rather  than 
jt-U-ycj  V  '      C^11  a  continuous  three-dimensional  model,  the  phase  change 
takes  place  on  an  interface  on  which  U  -    V      >a    ~      '  £     .  The  option 
here  arises  from  the  fact  that  in  the  one -dimensional  model,  the  depend- 
ent variables  are  averages  over  portions  of  stream  cross -sections .)  This 
has  worked  out  better  than  the  alternative  because  the  mass  -£-  A  t       is 
normally  taken  to  be  due  to  some  part  of  the  energy  transferred  into  the 
liquid  in  time  At.  Hence,  in  a  small  finite  interval  At,  there  is  less 
danger  of  removing  more  energy  from  a  stream  than  is  available.  Thus, 


S--h„(tt*.+  1,',)/A 


(ik) 


A  considerable  effort  was  expended  in  integrating  this  system  by 

[3] 
the  Lax  method   (Appendix)  although  it  fails  to  fit  the  model  of 

"Conservation  Form"  by  the  presence  of  the  term  -p  (Az),   in  (k)   and 

the  analogous  term  in  (5).  The  instability  which  plagued  this  approach 

was  probably  not  due  to  this  fact,  however,  but  was  due  to  the  fact  the 

system  may  become  non-hyperbolic,  i.e.,  it  may  have  complex -valued 

signal  speeds,  as  is  shown  below.   Since  the  Lax  finite  difference  form 

is  subject  to  the  "domain  of  dependence"  restrictions  associated  with 

the  maximum  signal  speed  in  the  medium,  a  calculation  of  signal  speeds 

was  made.  By  expanding  the  (2V2T)  system  to  a  quasi-linear  system  of 

the  form 

a .  ■    u  +    6,     a.       =0 

\j       tit  <</     h* 

[2] 
the  signal  speeds  X  are  found  as  solutions  of  the  system: 


Following  is  a  6th  order  algebraic  equation  in  \  whose  coefficients  are 
functions  of  the  flow  variables  (but  not  their  derivatives): 

o  =  (A  -  u)  (\  -  v)\(X-of[^  -  v  -  (A  -  v)3] 

-^(x-v)[^v  +  »Qi-  v)]  (15) 


where 


CO  = 


S  2 


*>*< 


dr 


(cO  >  c) 


(~l   >  o) 

(v  <  o) 


(16) 


In  general,  £-<->  will  be  small,  since  F>  .OS  will  be  taken,  and  will 
vanish  in  the  case  of  all  gas  where  7  -  / .  When  to  is  small,  and 
O  i=-   V  }   we  find  the  roots 


A  = 


where 


94. 


\ 


V 

u 


Uf  U)  v*    <f>     +       O^) 


(U-  v)((~.-v)V  +  *  °) 


'/? 


(17) 


(18) 
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Then  as   U  ->  V 


4 


,^-c-  v 


v  (u-  v) 


which  is  imaginary  for 


V(u-  v) 


<  o 


For  large  <-0  , 


\    -   < 


o 

V 

V  (i-  "Vj 


(19) 


V   _    y-  ^  K/  all  real# 

On  the  other  hand,  if  there  is  sufficient  friction  between  the  streams 
so  that  | V  -  U|  is  not  large  at  interior  points  of  the  (x,t)  domain  (as 
must  be  expected  for  small  <-o  ),   the  above  limits  do  not  apply.  Let  us 
assume  the  friction  is  large  enough  that 


U  ~    V    +       * 


CO 


(20) 


where  0~     is  bounded  as  <-*}  -*  °  .  We  obtain  the  roots : 
f 


V 

u 


A  - 


(21) 
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u 


+    nJ 


7Ul-») 


(i±\f^7 


sU-~  \J 


&(*$ 


which  are  all  real  for  small  co  at  least .  No  attempt  was  made  to 
determine  how  stringent  is  the  restriction  on  the  interface  friction 
and  initial  values  of  the  functions  imposed  by  the  requirement  of  real 
roots,  and  a  study  of  this  (2V2T) -system  was  postponed  in  favor  of  the 


(1V2T) -system.   (improperly  posed  initial  value  problems,  as  would  "be 
one  with  complex  roots,  occur  in  unstable  problems,  as  for  example  the 
flow  following  the  collision  of  two  liquid  masses). 

3 .  The  Reduced  System  (1V2T) 
Henceforth,  we  shall  discuss  only  the  system  with  a  single  velocity 
V  obtained  by  adding  (3)  and  (k)   and  setting  (J=  V     throughout.  We 
obtain  an  equation  for  the  conservation  of  total  momentum  which  is 

(A{  V)    t  (A(V\  A4)r4Ajt  -  S,  (22) 


f 


5         I-    i      +       f    (l-   g)  (23) 


^3   "     A  f   f   -   TJ  4     -    *j  £j 


AQ- *)*.„- A  *~    -  S  *  S 


(sk) 


Tv~n*?i 


3(  3i 


The  resulting  system  (1V2T).,  is  entirely  in  "Conservation  Form"  (A, 
a  given  function:  see  Appendix)  and  was  programmed  in  the  lax  scheme. 
Difficulties  are  encountered  when  area  variation  occfurs  due  to  the 
structure  of  the  Lax  Derivatives .  An  alternate  (lV2T)_.  system,  also  in 
"Conservation  Form" is  obtained  by  removing  A(x)  from  the  time  derivative, 

Let       B(x)  =■     /A  QO;    then 

ZGi>-    A* /a    -  -B\/B  (25) 


zt  -((/-») y)  -o-*jv2  --  S,b/s 


(26) 


(fi)     +   (rzV)      +    nVZ   -■<$   B  (27) 


J  t  J  X 


10 


(28) 


(rZ(^V\E))t+(^V(iV',E,^)l  ^ 

\i-di  V r  F))t  +(b-i)V(j  v'*f+  i^)) 


(30) 


This  was  programmed  in  the  lax  scheme  (using  ^  =  ^-=  °)  •  The  result 
■was  entirely  unsatisfactory  due  to  the  essential  averaging  process  implicit 
in  the  system.   It  appeared  that  the  CIR  (Courant-Isaacson-Rees)    method 
would  avoid  much  of  this  spreading,  and  could  yet  handle  shocks  by  inclusion 
of  the  concentrated  pseudo -viscous  diffusion  {^i-c,     ^p)  •     These  terms 
must  he  considered  as  source  terms,  however,  and  not  a  part  of  the  signal- 
generating  structure  basic  to  the  CIR  methods .  The  presence  of  the  dif- 
fusion terms  required  the  inclusion  of  an  upper  bound  on  the  quantity 
^/fcj  *)    based  on  stability  in  parabolic  systems.     Only  when 
strong  shocks  are  present  is  this  term  comparable  to  the  maximum  signal 
speed  limit.  The  term  used  was 


f 


•  /    I ]     /shock  diffusion! 

'?!]'[32t*IV(x*Ax)-  V(x)\)  I   condltion   J 

A  t   s<  {  (31) 


t'l  M  \    M/ 


/domain  of  dependence! 
"|     condition 
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The  equations  treated  by  the  CIR  method  are  the  quasi -linear  forms 
derived  from  (lV2T)>_,  and  called  (lV2T)_: 

F    f  Vzx    -  ((-*)    1/   =  D  (32) 


**  *   ^.  *  f  Vx     -  D,  (33) 


Vt  *    1/1/     +    i(PEM+  Pr    )-   b 


(3*0 


(35) 


f 
jt  ■>  x  ^f  l>£      JX  ^     >x/ 

where  for  convenience  we  have  "broken  the  sources  into  the  following 
forms: 


(36) 


D.   -   5u   *     D.akffJ*    Pt3(efiJ 


(37) 
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D„  -  V2  0-*) 


D    «- 

b    -  i 

(3                0 

K-~ 

rV  7 

»„-- 

33                         0    2 

*■  X 

A 


3/ 


"  7"  KVj  f  T*0*  4 


/) 


33 


K^3 


&,-   £(?.",  -  ?,*,-^-<^*2 


*   7   (n  * ,  *  *,  x*  ■+  ^  J 


0+i  -  [z-.f-  i+)/r 


D 


+  3 


ra* 


4"        ff^^^-I.O-yr-^^ 


r 


V 


( 


+  Y  Cr***  *■  **«* 


A) 


D 


4l 


«        <rY/-*J 


0    - 


/ 


v. 


vT3 


J^/-*J 


(r- 


E  -  ^ 


(38) 
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where 


v- 


Bi. 
Sjl. 


(39) 


(ho) 


T 


*'"(**)'(&)'  (& 


2 


<  O 


>  o 


(M) 


"ft* 


afs- 


o 


<^m 


"2V- 

8>j 


(te) 


where  a  is  a  number  between  1  and  2  (in  the  case  of  air  at  standard 
conditions,  this  has  the  effect  of  spreading  steady  shocks  of  moderate 
strength  over  three  to  four  net  intervals) . 


14 


This  system  is  of  the  form 


«■  bLt    U* 


-   D. 


6  -  /,•-,  ^ 


(^3) 


(summed  on  repeated  index)  with 


fel>    - 


/- 


V 

o 

c'3 

0 

V 

c" 

o 

c" 

V 

/c*' 

c41 

CHS 

c 


fl 


V. 


C 


n 


c 


r? 


o 

\ 

o 

o 

o 

^-  3  + 

C 

o 

C*t 

o 

crt 

\/ 

j 

m 


r 


c'3  - 

-  0-*) 

c"- 

r/i 

C*a- 

-  v/ 

C"- 

A>  0 

C*'- 

-•% 

C*J  = 

^-/>ca/r/- 

;J 

cf/. 

-  \A£  /f»C»- 

■^ 

c"<  = 

-C^fjV 

,c" 

c"- 

M 

^r4 

/.       /"  1  I  / 

^-     s? 

f 

c 


Ct5) 


V. 
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We  shall  assume  the  perfect  gas  thermodynamics,  which  gives 


where  \     is  the  ratio  of  specific  heats .  The  signal  velocities  now 
satisfy  the  5th  order  system: 


f(A-  vT  - 


o 


7 


V- 


PJ  P* 


While  it  has  not  been  proved  that  the  cubic  has  three  real  roots  for 
all  realizable  values  of  the  flow  variables,  these  roots  were  computed 
from  the  explicit  form  for  cubics,  and  no  complex  roots  were  encountered, 
(the  program  stops  in  the  event  of  complex  roots).  By  Descartes'  Rule  of 
Signs,  it  is  readily  proved  that,  if  all  5  roots  are  real,  there  are  at 

least  k   roots  (X       ,  ...,  X^   )  having  the  same  sign  as  V,  and  at  most  1, 

(5) 
X       ,  with  the  same  sign  as  (-V) .   In  the  case  of  "supersonic"  flow,  all 

X    have  the  same  sign  as  V.  As  y — » <*>,  only  the  first  two  of  the  cubic 

terms  persist: 

(K-V)+-o         ,  ana    (A-vJ-  "&-?(<-*}  v 

The  cases  Z -* /  lead  to  the  pure  gas  forms, 

ACjj  =    V 

and  (48) 


C*J       i  <rj 


Aw    y    ~    y  ±  fc+M 
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while  the  near-liquid  case  H-->  o  gives 

A  ®  =    — c  WW    1/ 

and  (49) 

giving  infinite  signal  speeds  as  2  -=»  °  .  However,  we  shall  in  practice 
take  z  >  z  min  =  0 . 05 •  Then 


<  3  6    yj  (^+-  Ivl) 


giving  a  lower  signal  speed  than  pure  gas  alone  (since  *p-  >  20,  normally) 
Thus  with  a  reasonable  restriction  on  z,  z  >  0.05,  the  limiting  signal 
velocity  will  normally  be  that  in  the  gas . 

If  all  the  eigenvalues  (signal  speeds)   A   of 

l\(W    -      b1*  I  -    o  (50) 

are  real,  we  assume  there  exists  a  non-singular  transformation  satisfying 
(p) 


? 


1  (^'-x(e)  r*)-°  (5i) 


This  has  not  been  proved  rigorously,  but  the  following  argument  makes 
it  a  reasonable  assumption  for  this  physical  system.   If  the  h   x  k 
matrix  obtained  from  (58)  by  deleting  row  2  and  column  5  is  non-singular, 
then  of  course  the  5  x  5  matrix  (58)  will  also  be  non-singular.  This  is 
related  to  the  fact  that  equations  (32-35)  of  (1V2T)  do  not  depend  on  F. 
For  the  k   x  k   system,  however,  a  non-singular  2   *  exists  in  general  when 
the  eigenvalues  are  distinct,  which  is  "almost  always."  A  small 
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independent  change  in  V    could  get  one  out  of  the  special  case,  however, 
yet  cannot  he  expected  to  change  the  mathematical  structure  of  the  system. 
Equation  51  defines  a  set  of  left  eigenvectors  £    corresponding  to 
X   ,  (p  =  1, . .  .,5)  .    Then  applying  this  transformation  to  (28),  we  get 


or  (52) 


[3] 

This  is  the  so-called  canonical  form  for  hyperbolic  systems   3    in  which 
each  equation  (index  p)  has  derivatives  in  a  single  direction  in  the 
(x, t  )  plane  defined  by  the  "characteristic  lines"  /    : 

r  p 

where  V   is  a  parameter  along   / 
Then 

u  -■  r   \ 


<-J    .    +       A        <u 


t  >* 


so  that 


JVP 


zpt 


t 


J  f 


u  =2 

^p 


1  &       (p-i.-,<r)      ^ 


is,  for  each  p,  an  ordinary  differential  equation. 

It  is  in  the  canonical  form  (^3) (using  (V?))  that  the  CIE  finite 
difference  scheme  is  analyzed.  When  the  coefficients  and  initial  data 
satisfy  certain  smoothness  conditions,  then  as  A  x  -=»  o  f   the  finite  dif- 
ference scheme  is : 


z 


zpL\-^- 


A  t 


+  £zpi A®[±  u>Hx+A\t)~  ^(^))_  7  /y > 


(5*0 


id 


.       (  AJL1 

where  A  t    £  /v*^*\>    \    [\\   j    and.  one  uses  (  +  )  according  to 
Ik,  A  ]     l  - 


is  stable  and  converges  to  a  solution  of  (V3)  in 


some  finite,  Ax- independent  time  interval  AT  >  0  of  the  initial  line. 
While  theoretically  valuable,  the  analysis  does  not  give  useful  estimates 
of  AT  or  the  error  for  Ax.  Certainly  the  procedures  could  not  be  contin- 
ued beyond  the  formation  of  the  first  discontinuity,  and  might  fail  before 
that.  The  addition  of  the  non-linear  pseudo -viscous  terms  (Ul-^2)  has  in 
practice  permitted  continuation  beyond  the  formation  of  shocks,  but 
estimates  of  the  error  must  be  inferred  from  knowledge  of  the  correct 
solution  or  experimentation  with  different  values  of  Ax.  Because  of  the 
large  number  of  arithmetic  operations  involved,  a  minimum  Ax(AT)  exists 
for  a  given  digital  word  length  below  which  accumulated  round-off  errors 
exceed  the  truncation  error.  For  smaller  values  of  Ax,  all  significance 
in  the  result  will  be  lost.   In  the  problem  at  hand,  however,  it  appears 
that,  with  the  8-digit  word  available,  the  computing  time  involved  for 

small  Ax  is  a  more  stringent  restriction  by  far  than  the  round-off  error. 

IZ  ^  i  j 

Taking  Z7    to  be  the  inverse  of   2  *  , 

*''<  2''  -     *''''  Z*>     -  f'  (55) 

and  we  return  the  system  (52)  to  a  form  where  each  equation  gives  a 
single  time  derivative: 

u'    ♦  Z  2"V'Aw'i?''     -Dl  66) 

where  the  summation  is  explicitly  shown  because  of  the  multiple  occur- 
rence  of  the  index  p.  The  index  p  on   U     is  to  indicate  that,  in 

a  finite  difference  scheme  according  to  (5*0  >  one  should  use  right-  or 

\(P) 

left-hand  difference  approximations  according  as  A   is  negative  or  positive 
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Since 


i~P 


2  r/*"(V"-  Awf,)J?    - 

by   (51);   using    (55)  we  obtain 

ft""-  \vDzu-  o  . 


C 


(57) 


so  that 


Z 


*Q) 


is  a  right  eigenvector  corresponding  to  the  eigenvalue 


(*) 


\W,  s=   (1,...,5). 


(3) 


As   V-»  <?,  there  is  exactly  one  root,  Xu  ,  of  the  cubic  (V7) 
which  vanishes .  Write 

then  explicitly,  we  may  choose: 


/" 


/ 


o-*J 


-\ 


z^ 


•*-t/" 


?z\j(A^U3liiM0-?)]  -i(A[3)-i)     0 


r 


A  -  (A^  A(r)J 


^    a 


V 


(58) 


-^ 


corresponding  to  eigenvalues 


Av"  -    V.   I/. 


>>    ■> 


°V    A(4)    Aw 


(59) 
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The  first  approach  to  the  finite  difference  numerical  integration 
was  to  use  (58)  and  the  computed  inverse  Z.         in  the  finite  difference 
form  of  (56) .  Difficulties  were  encountered  due  to  widely  disparate 
orders  of  magnitude  and  "ill-conditioning"  in  (58) .  (The  computed  product 
generated  diagonal  elements  near  unity,  but  with  some  off -diagonal  ele- 
ments  of  the  order  10  .)  Rather  than  go  to  double  precision,  a  satisfac- 
tory alternative  was  derived  "based  on  a  fact  earlier  mentioned  so  that 


we  may  put 


£-Ouvo 


(Af»)  - 


S-eyY^/ 


(V)  (J-',  ■■;*) 


J 


We  may  then  rewrite  (56)  in  the  following  form: 


u>t  *■ 


-h  K 


t3 


(60) 


-Dc 


where 


oC   = 


P 


and   (       2 
to        x<5) 


U) 


I 
-/ 


Cv>  o) 

tV<c) 


(61) 


/     aU)  v  <  o) 

Or) 


(62) 


n  l 

21  )   are  left-and-right  eigenvectors  corresponding 

Specifically,  we  take 
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r 


H       -    V 


3+    ,.    4-1  11   ,  +*  .    (s) 


?'*~~  ^c^l  c"c~-  o"(c"  * 


S3 


_  (c^-x^x^-v) 


c 


"hi 


ss 

Z  =     O 


\ 


and 


I 


ir 


c 


43. 


C 


ti 


ir 


3S 


\{s>-  V 


c 


IS 


i  r 


13     s~  3f 


C"  C 


I 


s~sr 


i 


ir 


(63) 


(00 


C    C     1-   c     c 


V. 


c,2(^-v) 
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For  the  'subsonic  case  ,  Q   -   '  ,  X*  '  has  multiplicity  1,  so  2     and 
-7  i(r)  ^  •?/(*} 

are  orthogonal-   Orthonormality  gives  the  unknown  scale  jr     : 


zt 


z 


/(f) 


~  I       /     £     EW*I^,    I  (65) 


^-.   Auxiliary  Routines 

Primary  emphasis  so  far  has  been  on  constructing  a  stable,  control- 
lable and  flexible  scheme  "FLOV  for  handling  two -phase  flows  based  on 
the  hydrodynamic  model  (1V2T) .  Auxiliary  factors,  such  as  boundary 
conditions,  initial  conditions  and  sources  terms,  have  not  been  inten- 
sively studied  or  varied  except  in  so  far  as  they  have  directly  interacted 
with  our  primary  purpose.  To  the  extent  that  they  have  not  interacted, 
simple  representative  types  were  taken  for  the  auxiliary  factors,  and 
elaborated  only  to  a  degree  that  seemed  necessary. 

As  far  as  was  possible,  the  auxiliary  factors  have  been  isolated 
in  subroutines  to  improve  organization  and  facilitate  changes . 

In  addition  to  the  basic  program  (FLO^)  using  the  CIR  method  with 
equations  (60-64),  three  subroutines  are  used:  BNDRY,  INVC  and  KAPA. 

Subroutine  BNDRY  is  used  to  set  values  of  the  dependent  variables 
on  the  end  points,  so  that  only  interior  points  are  treated  by  FIDk. 
In  cases  I  and  II,  BNDRY  has  been  based  on  the  assumption  of  an  infinite 
reservoir  at  the  left  end,  with  source  free,  adiabatic  flow  from  the  .right 
end  at  atmospheric  pressure.   Some  studies  were  made  with  a  finite  tank 
BNDRY  to  determine  tank  pressure  -  surge  pressure  coupling,  but  this 
effect  was  determined  not  to  be  basic  to  the  second  surge  formation. 

Subroutine  HWC  computes  and  sets  the  initial  values  of  the  depend- 
ent variable  (including  A(x)  determined  from  a  functional  form)  at  both 
interior  and  end  points . 

The  subroutine  KAPA  computes  for  interior  points  the  source  strength 
jz^     and  ^o-  and  the  source  coefficients  U'  ■      of  (38)  .>  including  the 
friction  terms  (  T^   ~g       ),  the  heat  transfer  rates  (  *?  0  %,  «  %i)} 
the  wetted  perimeter  per  unit  area  (  Ke  K»    K ^     )  and  the  non-linear 
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pseudo -viscous  stresses  and  work  terms  (  /*0,      /y         /(]        )  via  (kO-k-2) . 

The  (  Ke     K  a      l<i  )   describe  the  character  of  the  flow,  K0    indicating 
the  degree  of  dispersion  of  the  media.   They  are  prescribed  functions  of 
the  flow  variables  and  geometry.   In  cases  (i  and  II),  we  have  taken 

K£   =  (j-z)    P  (M) 

where  P  is  the  pipe  perimeter, 

K     -  i    P  (67) 

9 


and 


k    -    3  «i 


(68) 


This  gives  volume -proportionate  wall  wetting  in  gas  and  liquid  streams. 
Numerical  experiments  with  annular  flow,  Kg  =■  P^      K-~  &}    K0  =  {j~2)    P      , 
showed  sharply  curtailed  At  as  a  —=>   O    because  of  large  heat  transfer 
into  small  quantities  of  liquid.   Flow  experiments  indicate  annular  flow 
is  not  the  common  mode  in  Gooldown,  but  that    °'^l  should  be  much 

larger . 

Tke  (  T a     T.   )   are  taken  as  normal  friction  stresses: 

r     =    r  Vlvl  f     /? 
7  (69) 

t,    -    s    V/v/  ft   /t 

where  J  -  t,  -        GZ3     for  Reynolds  number  vd.       *d,z       greater  than 
2000,  and  /   /,  «  -^- .  for  Reynolds  number  below  2000.   Here 


214- 


(70) 


£  " 


*    /w 


^-c, 


^   1/  fi  (x) 

For  the  equation  of  state  in  the  gas,  we  take 


f   =  A,   r   T  (t\r  o.3(lo)7) 


(71) 


A 


with  gas  and  liquid  energies, 


(72) 


F  -    Cvl  (Tt  -  Tj 

c   =-  o.  fi  (i°y 

We  take  the  following  linear  approximations  to  the  pressures  and  latent 
heat  from  data  for  N  ^     on  the  vapor  pressure  curve: 


and 


^a."  K,    '     /^^7"-  -4./3*(/0)'t    &ZIUQ*)'   T 


(73) 


and 


\,  2,  +£j  T  =  zj?(fo)f-  i^o(io)7  T 


Since  A  -  (£"-  FJ  ,   we  find     £*^  «  O.  Q  1 1  (f  o) 


(Ik) 


s,t. 


and    /  c  =  5  3  0  * 
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2 
We  take  the  heat  transfer  rates  per  cm  per  degree  Kelvin  from  wall 

*  *  * 

to  gas  Q      ,  from  wall  to  liquid  9  £    ,  and  from  liquid  to  gas  q         to 

^"  h  fix  1  '"* 

"be  a  simple  form  of  the  Dittus-Boelter  correlation: 
*       *  *   '    i ■  —  °-? 


1]-    %"  ^JBC*)   «. 


%]  -    %V  ^t  iBU)  & 


I 


(75) 


1*.  -    <"  d,i/86tl  ^  (76) 


and  put 


•/ 


A    ■=    o.i(io) 


a  (3^3 

Then  for  heat  sources  we  have 

f"j*7-    Vj'fc'" 

V 

5 

<i  %e   -     K*   1*  (t„- 

7-J 

U„?.  -   *.  %*  Cr*  - 

^ 

(77) 


Since  we  are  using  only  "bulk  heat  transfers  and  mean  temperatures, 
no  precise  estimates  of  evaporation  or  condensation  at  the  interface  can 
"be  made.  Such  an  estimate  would  depend  on  temperature  gradients  at  the 
interface,  and  the  relative  values  of  interface  temperature  and  interface 
saturation  temperature.  We  distinguish  three  sources  of  evaporation  and 
condensation:  (  jz,       ^c   )  at  interface,  (  ^a  ^j  )  at  wall,  (  -e^   •«•,  ) 
spontaneous  for  supersaturated  cases .  Finally, 


(*,,   ^>)    -  (  **>,;    **,  )      +     U^   ^aJ  *    (^    <**) 


(78) 
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Then  we  have  used  the  following  forms 


i 


O 


/ 


K 


0     *^0 


'-fc-zj 


{T<Tftt(f))     (79) 


^-^A  (Tj>T.tW) 


I 


Q 


-   T. 


(  z  *  Z«,J 


A     [Ts,t>rt 


(80) 


where 


%t  -  srv^c.  [  e-  Fj    o.Qrhof] 


(81) 


X. 


< 


o 


<2 


<t  ll(Tl-  77, J 


/I 


ft;*  77,  J   (82) 
fc  >  7;,, 


and 
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r 


o 


^u 


v. 


V~    K3  ^g  l_Jjl  '^'A 


(83) 


The  0~)    are  "efficiencies"  of  conversion  of  transferred  heat  to  evap- 
oration and  condensation.  Also, 

r 


X3  = 


o 


(*  >  2t?Ax) 


r<  *  Z.t 


r,o-*)h-vjA      b:>r,j 


m 


and 


^r 


*i*(r„r-TjA 


(85) 


We  have  used   CT  =  0,  \T  f  ^  -  f ,  ^    *tj    S~)       ,  and  <r  •=  0.  00/'  for 

(^  =■  3j  6).  This  latter  was  chosen  simply  because  it  gave  quantities  for 

(  ^  /cz      )  of  the  same  order  as  had  been  obtained  in  typical  cases 


for  (  -£a  ^ca  ). 


2  V 


In  the  calculation  of  (4l,  ^2)  the  calculation  of  <=r-"   was  modified 
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to 


f 


-  I. 


the  minuend  being  chosen  to  remove  the  part  of  the  derivative  due  to 
variable  geometry,  as  computed  from  the  steady  adiabatic  flow  equations 

(A(i-*)vlx-a 

(An- v)  ,-  o 


J   X 


In  the  equation  (7)  for  wall  heat  transfer,  the  heat  transfer  rate 
was  taken  to  be 


<A 


^■-  -77^^  =  ""^  (87) 


where  P„   is  the  perimeter  of  the  uniform  pipe,  S     its  thictaess,  ^ 

its  mass  density,  and  dpu_,its  heat  capacity. 

5 .   Controls  and  Blocks 
This  finite  difference  scheme,  called  FIDk,    is  finally  then  in  the 
form 

ml(*;tfAt)=  <^l(\t)  *  At(Dl(x,t)-  E'fx,  t)J  (88) 

At  least  three  complete  passes  through  the  net  of  points  in  x  are  involved 
in  FL04  for  a  single  time  step.  On  the  first  pass,  boundary  values  are  set, 
and  (At)   ,  satisfying  (3l)>  is  obtained.  We  also  require  At  to  be  less 
than  twice  the  value  of  At  used  on  the  previous  time  line .  D  and  E  are 
also  computed.  On  the  second  pass,  (u'J  s.    u  '  ( Xj  t  f  A  t)   is  found  at  each 


29 


x,    and  tested  to  be  certain  that  the  following  limits  are  satisfied: 
(  l2r-2l     «    (OLl)-   1*1 

\r+  -  r /  t   (GLR)-lrl 

[\J+-y\   ,<   ^w((uuJy(uLuJ.|v/J  (89) 

W-Tj  <  (ullJ-/t-J 

Tested  for  positivity  is  the  argument  of  the  logarithm  in  (92).   If 
any  of  these  tests  fail,  we  reset  At  =  2"At,  and  repeat  the  second  pass . 
In  the  event  that  the  computed  value  z  falls  outside  the  range 
7    x    2:     $    z? _  .  .  ,  we  set  z+  to  the  nearer  of  these  two  limits,  and 
use  equation  (32)  in  reverse  to  solve  for  the  --o   or  -£.  (positive- 
valued)  required.  This  latter  value  is  used  in  the  remaining  equations. 
In  case  I  and  II,  we  have  used 

U  L2  -    OL&-    UUO  -  (jLt-  ULL--  G.  I 

ULV=     <3.  /  Uo) 

-to  (90) 

In  the  third  pass,  the  new  values  of  the  dependent  variable 
(z+(t,x),  etc.)  are  stored  as  old  values  (z(t,x),  etc.)  in  preparation 
for  computing  a  new  line.  BWDRY  is  used  here. 


In  view  of  the  fact  that  1  -  z  must  be  permitted  to  become  very 
we 
replaced 


-  <r 
small,  we  have  used  /—  £  =  C     in  all  algebraic  expressions,  and 


A  2  -    F       At  (91) 


by  its  exact  equivalent 

A  c^=  -J^  ( I  -  cf  At  ■  F\ 


(92) 


and  similarly  for  x  -  derivatives 
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6.  Computation  Summary 

Computations  were  carried  out  on  an  IBM  7090  with  an  8-digit  floating 
point  mantissa.  The  finite  difference  scheme  "based  on  (60-6k)   appears  to  be 
stable  and  to  give  solutions  of  the  proper  character  during  the  growth  and 
early  stages  of  decay  of  the  first  surge,  a  real  time  of  about  2  seconds . 
In  Case  I,  a  peak  pressure  of  17  atmospheres  was  calculated.  For  these 
computations,  Ax  =  290  cm  was  taken  (20  stations  in  the  pipe  length), 
yielding  a  computing  rate  of  0.6^  flow  seconds  per  minute  of  computing. 
Since  flow  times  of  the  order  of  10  seconds  or  more  are  desired  (l6  minutes 
computing),  it  does  not  seem  reasonable  to  reduce  Ax  significantly,  espe- 
cially since  cutting  Ax  by  two  quadruples  the  computing  time  for  a  given 
real  time  interval. 

After  the  first  surge,  liquid  was  flowing  back  into  the  tank,  and  gas 
was  flowing  rapidly  from  the  down  stream  end  of  the  pipe.  The  pressure 
fell  rapidly  until  it  began  to  near  the  tank  pressure.  At  this  level, 
the  pressure  no  longer  fell,  but  was  sustained  partly  by  boiling  of  a 
non-diminishing  liquid  residue  in  the  pipe  and  by  non-diminishing  gas 
energy  and  gas  density.  A  quasi-steady  state  thus  resulted,  and  no  second 
surge  evidenced  itself. 

The  reason  that  energy,  pressure,  and  liquid  phase  are  maintained 
in  spite  of  the  rapid  transport  of  these  quantities  out  of  the  pipe  at 
each  end  can  be  seen  by  writing  the  equations (60)  in  the  form 


U>t       * 


^l^  As-) 

=  -  x(3  Z        A 


(93) 
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and 

P L   ^     J  (Case  I)   (9^a) 


where  we  have  approximated  the  difference  in  right-  and  left-hand  deriv- 
atives by  a  second  derivative  within  the  accuracy  of  the  finite  difference 

approach.  The  equations  (93)  a**e  diffusion  equations, 

v  tr)  ^iW  is)  i 
-o<"AZ_    £     >  O         ,  the  diffusion  terms  in  the  right-hand  side 

tending  to  smooth  out  irregularities  in  the  functions  and  their  derivatives 

Under  the  smoothness  conditions  of  the  CIR  proof,  these  terms  vanish  as 

/S.  X  — >  O  ,  but  remain  large  enough  to  stabilize  the  numerical  method. 

In  the  computation  described  above,  however,  these  terms  dominate  the 

convective  terms  on  the  left-hand  side  when  the  pressure  has  fallen  to 

nearly  tank  pressure.   This  diffusion  also  probably  accounts  in  part  for 

the  earlier  surge  peak  times  encountered.  This  numerical  diffusion  is, 

of  course,  a  physical  interpretation  of  the  effect  of  truncation  involved 

in  the  finite  difference  equations . 

A  few  attempts  have  been  made  to  circumvent  these  difficulties,  and 

allow  for  the  formation  of  multiple  surges  without  reducing  the  stability 

of  the  system.  The  only  moderately  successful  one  was  to  introduce 

multipliers  which  tended  to  reduce  the  diffusion  when  it  was  already 

small : 


/fu:.).-(«U( 

/  c     \    I       [f   i    \    I  (Case  I3:)  (9^) 


Fl  = 


Ft  i 

*  /    ,  and  t~      -  I      if  the  right-  and  left-hand 

derivatives  have  opposite  signs.  Using  this  form  (Case  II)  a  small 

second  surge  was  observed,  while  the  peak  pressures  in  the  first  surge 

reached  only  11  atmospheres . 

Because  so  severe  a  flow  change  appears  to  result  from  adjustment  of 

purely  computing  parameters,  it  might  on  first  sight  seem  to  invalidate 
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"both  results,  "but  the  diffusion  is  unquestionably  overvalued  in  each  case. 
Any  scheme  which  retains  stability  while  reducing  the  diffusion  should 
give  improved  results .   It  might  well  be  worth  testing  other  ideas  in 
this  regard.  Two  alternatives  are  (l)  greatly  reducing  Ax,  giving  a 
practical  computing  method  only  for  problems  in  a  very  small  geometry 
or  a  very  short  time,  or  (2)  attempting  to  program  a  Iagrangian  scheme, 
in  which  the  diffusion  required  for  stability  might  be  much  less,  and 
better  liquid  control  might  be  maintained.   It  does  seem  that  because 
of  the  high  mass,  energy  and  momentum  density  of  the  liquid,  a  good 
representation  of  liquid  movement  is  essential  to  a  meaningful  computing 
scheme . 

It  should  he  noted  that  in  (93 )>  *,   the  factor  in  the  liquid 
continuity  (z)  equation,  can  be  set  to  zero,  or  even  made  negative  with- 
out interfering  with  the  stability  of  the  scheme,  so  long  as   o  <  2-  <  I 
is  retained.   In  some  tests  here,  further  improvements  have  been  made  in 
avoiding  liquid  diffusion,  but  results  are  not  conclusive.  Making  F 
dependent  on  flow-direction  and   o~x~'  >   has  also  been  tried  inconclusively. 
Some  realism  might  be  gained  here,  in  reducing  the  requirement  of  conti- 
nuity in  the  liquid,  and  emphasizing  transport.   (indeed,  for  small  z 
with  droplets,  z  -  waves  cannot  propagate.) 

Cases  I  and  II  were  computed  with  the  initial  values:  (CGS  system) 
Ax   <x  $  J~  A  x  J~A  x  <   x    *  Jo  A  X 

V  o  a 

c.S'oGs-(io)  o.  i  o  i  J  (l  o)1 

a.  «W  /  H  (t  o)~'  a.  (i  3.  &  ( i  o)~* 

c  ot~  o.  ??  ??r 
7c*  K.  3  oc  °  K. 

T£  70°K.  t,fc  K. 

X  7  e '  K.  3oc*  K. 
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with  the  single  addition  /u  (  4  A  X^  oJ=  ~J O    K.   to  admit  heat  leakage. 

Another  variation  was  made  on  FI04  (using  F  =  l)  to  test  a  model  of 
valve  opening  time  T  and  to  attempt  to  match  computationally  the  strong 
dependence  of  peak  pressure  on  valve-opening  time  observed  experimentally. 
While  not  thoroughly  understood,  this  dependence  may  "be  due  to  liquid  jet 
formation  at  the  valve .   If  this  were  so,  the  effect  would  he  due  then  to 
spatial  variation  which  cannot  be  expected  to  come  out  of  a  one -dimensional 
analysis .  We  used 

B(*,t).   BM  =p       C°  «  *  <  tJ 

B(x,t)'~  S(x)  (t  ±T)  (95) 

An  increase  of  peak  pressure  did  occur  as  opening  time  increased, 
hut  this  trend  did  not  reverse  itself  for  large  opening  times  as  is 
observed  experimentally.  Again  the  liquid  diffusion  dominated  transport 
when  the  valve  opened  slowly;  the  liquid  that  was  thus  introduced  down- 
stream of  the  valve  boiled,  and  the  valve  restrictions  prevented  the 
required  back  flow  which  would  have  prevented  the  sky-rocketing  pressure 
at  10  seconds.  Furthermore,  the  valve  was  represented  by  a  transient 
area  variation  which  was  not  small  and  stationary  as  was  assumed  in  the 
derivation  of  the  equations . 

Again,  finer  difference  nets  would  have  permitted  better  valve 
representation  without  the  large  area  variations  introduced  here: 
Valve  Opening  Time  (Sec)     Peak  Pressure  (Atmosphere)     Time 

o.i  13.8  0.30 

1  15.9  0.U9 

3  19.3  0.59 

5  39.0  0.43 

10  ? 

(See  figures  3>  ^>    and  5«) 

In  Figures  1  through  5 >   "the  pressure  -ft     ,  volumetric  quality  z,  and 
velocity  V  are  plotted  at  each  of  the  twenty  pipe  stations  for  times  given 
by  the  baselines  (arrows  on  time  scale) .  Above  each  baseline,  h      is 
measured  in  atmospheres  (absolute )/o.l  inch,  z  in  units/inch,  V  in  1000 


3^ 


cm/sec/O.l  inch. 

Figure  6  shows  the  pressure  in  atmospheres  for  Case  I  (lower)  and 
Case  II  (upper)  as  a  function  of  time  for  station  7  (x  =  7  A  x),  the 
location  of  the  pressure  peak  in  each  case.  The  circles  are  at  the  points 
in  time  for  which  results  were  printed,  less  frequent  in  Case  I. 

Figure  7  is  an  experimental  curve  given  in  a  private  communication. 
Note  the  spread  in  experimental  data. 

Clearly,  the  worst  defects  of  the  FLOU  model  are  in  time  of  first 
peak  and  peak  decay  rates  .  While  perhaps  these  could  be  improved  "by 
modifying  the  geometry  B(x)  and  improving  the  empirical  heat  and  momentum 
transfer  formulae,  it  is  certain  that  the  numerical  diffusion  is  the  major 
culprit.  With  the  present  costs  of  high  speed  computing,  it  is  not  prac- 
tical to  use  Ax  small  enough  to  evaluate  the  model's  potential  effectiveness 
in  predicting  the  Cooldown  process. 
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7.   Glossary 
Physical  quantities  are  measured  in  CGS  units.   Subscripts  g,£,o 
refer  respectively  to  gas  (or  wall-to-gas),  liquid  (or  vail -to-liquid) 
and  interface  (or  liquid-to-gas),  and  are  replaced  "below  by  the  generic 
subscript   V  .  Units  are  shown  in  square  brackets . 
a.  .:  coefficient  matrix  of  time  derivatives  for  quasilinear  hyperbolic 

system . 
A  :  gas  constant  in  state  equation  for  perfect  gas.  (eq.7l) 
A(x) :  cross -sectional  area  of  pipe  (A(x,t)  used  in  study  of  effect  of 

valve -opening  time) .  [cm  ] 
b   :  coefficient  matrix  of  x-derivatives  for  quasilinear  hyperbolic  system. 
B(x):  l/A(x).  [cm-2] 
/c  }  ^o  . :  mass  rates  of  condensation  per  unit  length  of  pipe,  (eq.78) 


cm.  fee. 


[orgs     ] 

C   :  specific  heat  at  constant  volume.  ['    ~  -1  ^  j 

I  e^l^   1 
C   :  specific  heat  of  pipe  wall.  I    v  c^' 

D.,  D. . :  source  terms  for  (1V2T)  .  with  partition,  (eqs.37>38) 

-C    ,-C>.:   mass  rates  of  evaporation  per  unit  length  of  pipe,  (eq.78)  I'^Zri^'c  I 

E:  specific  energy  of  gas.  (eq.72)  /   '  *j  J 

1 
F:  specific  energy  of  liquid,  (eq.72)    I *—~  j 

F  :  a  non-linear  diffusion  coefficient,  (eqs  .9^-a,9^b)  [l] 

F*:  typical  time  derivative,  (eqs. 91*92  only) 

f   :  friction  factor,  [l] 

Q  (x) :  component  of  gravitational  acceleration  along  pipe .  /  - — a  I 

/ON     l^JJL-     1 

h   :  wall  heat  parameter,  (eq.07)   /  — t~    / 

WX  r  .,    <-~i<3  -* 


h^,  :  heat  conductivity.  [  ~  ;t 
k^:  i,B.  (eq.39)  [**»''] 


rt.Je.c-.   <\.  J 
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Jt^  '   perimeter  of  contact  on  pipe  cross-section,  (eqs.10-13)  [cm] 

/jo.  :   pressure.  J-£_r  J 

/p.     (t) :  vapor  pressure,  (eq.73)  I     I J 

P(r,E):  general  equation  of  state  for  vapor,  (eq.8) 

P  (x) ,    r     :   pipe  perimeter,  and  its  value  on  uniform  section,  [cm] 

J  /  e^x     ] 

q*  :  heat  transfer  rate  per  area  per  degree,  (eqs  -7^,76)1  Cma  jcc   o^-  / 

q*-*:  constants  in  Dittus-Boelter  correlation.  (eq.75,76)  [l] 

q  :  heat  transfer  rate  per  unit  area,  (eq.77)  I ffr~: — T  J 

^n     :  pseudo -viscous  stress.  (eqs  Al,42)   /  —+z — 3 
r:  gas  density.  [  ~T 3  / 
ftcL   :  Reynolds  Number,  (eq.70) 
s:  liquid  density.   /  CfA* J 
S.  :  sources,  (eqs. 9-1*0 
Jy   . :  sources.  (eqAo) 


1 
t:  time,  [sec] 

T:  valve -opening  time,  (eq.95)  [sec] 

T  ,  (fa):   saturation  temperature,  inverse  of  /h       (T) .  [ °K] 

T  :  fluid  temperature .  [ °K] 

T  :  wall  temperature.  [ °K] 
w 

T  :  "base  temperature  for  eqs .  72.  [  °K] 

u  :  dependent  variable  vector:  (z,r,V, E, F) 

(u  ,    )_,  (u  ,  )_:  right-  and  left-hand  x-derivatives . 
x  'x  R'  v  'x  L 

few.   / 
U:  liquid  velocity.  [  "Xe  0  J 

ULZ,  UIR,  ULU,  ULV,  ULG,  ULL:  limits,  (eq.89) 

V:   gas  velocity.     [fzZ.i 
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x:  coordinate  along  pipe  axis,  [cm] 

y:  quality  (mass  of  gas  per  mass  of  fluid  mixture),  [l] 

z:  fraction  of  pipe  cross-section  containing  gas.  [l] 

Z:  pipe  expansion  rate,  (eq.25)  [cm   ] 

ZMAX'  ZMHT:  lDOUnds  on  admissible  z.  (eq.  90) 

z  ,  Z   :  matrices  of  left  and  right  eigenvectors.   (eqs  .51, 57>63, ^4)' 

Y   :  ratio  of  specific  heats  in  gas  phase,  [l] 

d  :   pipe  wall  thickness,  [cm] 

<3   :  Kronecker  delta . 
Ax:  fixed  space  increment,  [cm] 
At:  variable  time  increment.  (eq.3l)[sec] 
A  :  fluid  signal  speeds .  /  J^  J 

\*:   latent  heat  of  vaporization.  (eq.74)  I  ~  ?  '  / 
^r  :  viscosity.   I'ZZiZc.  J 
^  :  (eqs. 16,46)  ["^.J 

v  :  (eqs. 16,46)  /e-^j 

p  :  mean  density,  (eq.  3l)  /  "^ — i  / 

p   :  pipe  wall  density.  /-5 — S  / 
\   V  I    C  m  .  J 

sr  =  -  £^  (i-z)  [1] 

CT". :  coefficients  in  formulae  (79-83).  [l] 

T     :    shear  stresses,  (eq.69)   /   ^vT*~  / 

to  :  (eq.l6)  [l] 

(2V2T) :  system  of  conservation  equations  (l-8)  allowing  each  fluid  a 

velocity  and  temperature . 

(1V2T). :  system  of  conservation  equations  with  one  velocity  obtained  by 
replacing  (3,4)  in  (2V2T)  by  (23,24) 
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(1V2T)  :  system  of  modified  conservation  equations  (26-30)  with  one 
B 

velocity. 

(lV2T)p:  system  of  quasilinear  equations  (32-36)  with  one  velocity  used  in 
FL(A  numerical  evaluations . 
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Appendix 
The  Finite  Difference  Method  of  lax. 
Lax's  method  assumes  a  system  of  hyperbolic  partial  differential 
equations  in  "Conservation  law"  form: 


where,  for  our  cases,  (r y  <?  J    h    )  are  algebraic  functions  of  the  flow 
variables  (  2-   f"j    V ^    T ~      7~^  ) .  The  computation  is  advanced  in  time 
from  a  given  line  of  values  given  at  equally  spaced  (2  A  x)  values  of  x 
on  a  staggered  net  by  substituting  in  (A-l)  the  following  expressions: 


A    ,      f\   ^    -fl(*J  iAt)-i[fl(x+A*A^fL(x-A*J)] 


(A-2) 


ji  At 

and 


gl"(x,tj-   ?L'f"*x,^-  f  (*-**-  ti 


(A-3) 


2  A  X 


After  computing  a  line,  the  flow  variables  are  computed  from  the 
f  (t  +  At,x) .  On  theoretical  grounds,  we  must  have  that  At/   be  less 
than  the  maximum  signal  speed  occurring  in  the  flow  at  time  t.  The 
averaging  resulting  from  the  average  used  in  (A-2)  is  necessary  to  the 
stability  of  the  method,  but  causes  too  great  a  diffusion  of  flow 
quantities  in  the  Cooldown  problem  to  be  useful.  There  does  not  appear 
to  be  a  rigorous  theory  of  convergence  and  stability  for  the  lax  method 
comparable  to  that  given  for  the  CIR  method  (39) • 
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Figure    I 
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FLOW     STRUCTURE     IN     X    AT 
A    SEQUENCE     OF    TIMES 

CASE    II:    F'   *l    (SEE    Eq.  94b) 
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Figure   2 
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1.2 


INFLUENCE   OF    VALVE    OPENING 
TIME    ON    FLOW    STRUCTURE 

F1    =  I,    T  =  O.I    sec     (SEE    Eqs.    94a,   95) 


Figure   3 
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INFLUENCE    OF    VALVE    OPENING 
TIME    ON    FLOW    STRUCTURE 

F1    =  I,     T=  1.0   sec.    (SEE    Eqs.   94a,  95) 
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Figure   4 
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INFLUENCE    OF   VALVE    OPENING 
TIME    ON     FLOW    STRUCTURE 

F1  ■  I,    T=3.0   sec    (SEE    Eqs.  94a,  95) 
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Figure    5 


PRESSURE    VERSUS   TIME  AT    A 
POINT  20    FEET    FROM   THE    VALVE 


case  n      F'gfe  i 

(SEE    Eq.    94b) 
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CASE    I        F1  =  I 

(SEE    Eq.  94Q) 


Figure  6 
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TYPICAL    EXPERIMENTAL    CURVES  OF   PRESSURE 
AND    LIQUID    VOLUME    FLOW-RATE    DURING    FIRST   SURGE 
CYCLE     FOR   ONE    DRIVING    PRESSURE 
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fluid:   liquid  nitrogen 

TRANSFER    LINE:    200    ft.  X  0.75  in. 

0.  D.  x   0.67  in    I.D.    COPPER, 

VACUUM     INSULATED. 

INLET   DRIVING    PRESSURE    62.5   psig, 

TEMPERATURE     75.7°   K, 

VALVE     OPENING    TIME:    APPROX.    1/2    Sec 

EXPERIMENTAL 

COMPUTED 


DRIVING     PRESSURE 


(62.5    psig) 


Figure   7 
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